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BCS superconductivity is explained by a simple Hamiltonian describing an attractive pairing in- 
teraction between pairs of electrons. The Hamiltonian may be treated using a mean field method, 
which is adequate to study equilibrium properties and a variety of non-equilibrium effects. Nev- 
ertheless, in certain non-equilibrium situations, even in a macroscopic, rather than a microscopic, 
superconductor, the application of mean field may not be valid. In such cases, one may resort to 
the full solution of the Hamiltonian, as given by Richardson in the 60's. The relevance of Richard- 
son's solution to macroscopic non-equilibrium superconductors was pointed out recently based on 
the existence of quantum instabilities out of equilibrium. It is then of interest to obtain analytical 
expressions for expectation values between exact eigenvalues of the pairing Hamiltonian within the 
Richardson approach for macroscopic systems. We undertake this task in the current paper. It 
should be noted that Richardson's approach yields the full set of eigenvalues of the Hamiltonian, 
while BCS theory yields only a subset. The results obtained here, then, generalize the familiar BCS 
expressions for, e.g., the electron occupation or pairing correlations, to cases where the spectrum 
of excitations diverges from BCS theory, for example in cases where the spectrum exhibits multiple 
gaps. 
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I. INTRODUCTION AND RESULTS 



Our basic understanding of superconductivity is informed by the mean field solution of the pairing Hamiltonian, Eq. 
(Ill), given by Bardeen-Cooper-Schriefferi^ (BCS) some half a century ago. The mean field solution does extremely 
well in many respects because of the fact that, in essence, the condensate interacts equally strong with all energy 
levels that participate in superconductivity. Thus a macroscopic system is well described by mean field theory. There 
are some caveats, though. The BCS expression for the eigen-states does not capture all possible eigen-states of the 
pairing Hamiltonian. For example the BCS expression predicts a single gap in the excitation spectrum, while as it 
turns out from an exact solution, any number of gaps may appear in the spectrum. This is known because Richardson^ 
has solved the pairing Hamiltonian exactly. These unexpected eigen-states, which correspond to an unusual spectra 
of excitations do not play an important role in equilibrium - a stroke of luck for BCS theory. Nevertheless they 
may appear in out of equilibrium situations, when the system is far enough from equilibrium^, this is due to certain 
quantum instabilities encountered in non-equilibrium superconductivity^. 

In this paper we study more closely the consequences of these multi-gapped eigen-states. We make use of Richard- 
son's exact solution and Slavnov's formulai^ as applied to this model^ in order to compute correlation functions of 
the pairing amplitude and the occupation number at different energy levels for a macroscopic system analytically. 
It should be noted that Richardson already derived expressions for correlation functions, which are specific cases of 
Slavnov's formulas^ deriving their thermodynamic limit for the case of one gap and for the expectation value of a single 
operator?. Our a result is a generalization of that work, giving expressions for any number of operators both in the 
cases of one or two spectral gaps. Similar approaches were used in [lO and fill using numerics and focusing rather on 
mesoscopic systems. The relevance of Richardson's exact solution to the pairing problem in mesoscopic systems was 
pioneered in Ref. (l^ . Our results, dealing with a macroscopic system, agree with BCS theory, when there is only one 
spectral gap. Indeed, the expectation value in that case is given by Eq. below, (where i?2(0 = ~ + ^^); 
an expression which is familiar from BCS theory. 



II. THE RICHARDSON SOLUTION 



BCS superconductivity is captured by a model Hamiltonian, designed to include only those features, that are crucial 
to the existence of superconductivity. The Hamiltonian includes a free, bilinear, part composed of L single particle 
levels and an interaction part which scatters Cooper pairs: 

(Jj ^± l<l <L 

Here (j, and (j, — ) denote the quantum numbers of time reversed pairs. For example, if (j, +) denotes a state with 
wave number k and spin up, then (j, — ) denotes a state with wave number —k and spin up. We assume, for simplicity, 
that each level j is only doubly degenerate, where a indexes the two degenerate states, a taking + and — as values. 
Furthermore, we assume uniform level spacing Sj — £j-i — 5. 

It turns out that the model Hamiltonian is exactly solvable, namely, the eigenvalues and eigenstates can be found 
exactly. This was done by Richardson in the 60's'^, motivated by the Hamiltonian's importance in nuclear physics. 
The solution of the Hamiltonian is not trivial, and in fact the Richardson solution may be understood as a Bethe 
ansatz solution. Indeed, the Hamiltonian contains a non-trivial interaction which scatters Cooper pairs (namely 
time reverse pairs of electrons). Note however, that those electrons which singly occupy levels do not scatter. The 
levels which contain singly occupied states are then called blocked levels^ since no pair can scatter into them. The 
set of levels which are occupied by single electrons, together with the corresponding spins of the electrons, are good 
quantum numbers. A given set of such quantum number is termed seniority. For each given seniority, denoted by 
(ji, (Ti,j2, 172, ... jjM, ctm), one may define a 'vacuum' state \{ji:Cri]iLil i which contains no pairs, but only singly 
occupied levels: 

iu,^.}f£i)=n4.-jo) (2) 

i 

Another good quantum number is the number of rapidities, P. The total number of electrons in the state is then 
2P + M, since each rapidity contributes one pair of electrons, due to Eq. ([!]). There is a one to one correspondence 
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of eigenstates with given P and given seniority and solutions, {E^}f^^^, of the Richardson equations: 

Here U is the set of aU unblocked levels, namely V = {i\$k,i = j^}. The eigenstate is then denoted by 
\{Ei}iLi,{ji^'^i}iii) and is given explicitly by: 

a=l i^U 

The eigenvalue for the state is = 2E^ + J2j(fu ^J- 

We shall term the E^^s as rapidities characterizing the state K-B^jfLii {jii o'ilfii)- We shall term a state of the 
form @ a 'Richardson state' even if the rapidities do not satisfy ([3]). Namely, a Richardson state is an eigenstate if 
and only if the rapidities satisfy ^ . 



III. ELECTROSTATIC ANALOGY 



In what follows we shall be interested in computing expectation values in the Richardson model in the continuum 
limit. To do so, Richardson's equations ([3]) must be solved in different circumstances. The solution of the Richardson 
equations are facilitated by the fact that these have a convenient 2D electrostatic interpretation. The electric field at 
point w of a charge placed in at point z in the complex plane is given by — iEy = ^^^^ , which means that ([3]) may 
be interpreted as the condition of electrostatic equilibrium of the charges E^, which are assigned a charge +1, given 
the the position of charges of magnitude —1/2 at Si (only for i e [/ ) and a constant field, i, pointing in the negative 
X direction. 

Given this electrostatic analogy it is natural also to define the field h{z) = 5 {Ex{z) + iEy{z)) (here 5 is the level 
spacing). Explicitly, h{z) is given by: 



E 



1 



2 ^ £ - £, ^ f-E,, 



1 

g' 



(5) 



Here 5 = y . The continuum limit is taken by letting 5 — > while letting h{z) tend to a constant, except for certain 
arcs and a segment of the real axis. In fact the solutions of Richardson's equations have the following form in the 
continuum limit: as can be easily understood, there is always an electrostatic equilibrium position between any two 
adjacent unblocked Ej's. As the level spacing is decreased (5 — 0, one may create any line density of charges on the 
real axis, by placing rapidities on the real axis between adjacent unblocked levels or by blocking levels. This defines 
a coarse grained charge, p{e) on the real axis, given by 



p{e) = 5 ■ 



Ssie-E^)-]^ ~5s{e-ei) 



(6) 



here Ss is a function tending to a delta- function as 5 — > having a width much larger than 5, e.g. Ss{x) 



for some large number A. In addition to those rapidities, which are on the real axis, and lie between the e, there 
may exist truly complex rapidities. These, it turns out, arrange themselves on arcs. The arcs extend symmetrically 
around the real axis (because the Richardson equations are real, complex rapidities come in complex conjugate pairs). 
Suppose there are K arcs. We shall denote the two endpoints of arc j as iij ± iAj (here i is the imaginary unit). 

The field /i(^) consequently has a jump discontinuity on the real axis of magnitude p(e) and on the arcs, where it 
has some 0(1) jump discontinuity, which must be determined. Consider the endpoints of the arcs. Those rapidities 
that sit on the endpoints, must be in electrostatic equilibrium. A closer analysis shows that this is only possible if 
the field h{^) as ^ approaches the endpoint tends to 0. This is an intuitive result, since if h{^) would not approach 
the endpoints would feel a force that would move them. One concludes that h{£^) vanishes on the 2K endpoints of 

the arcs. Moreover, if we look at the average value of h{£,) across the arc (namely .Mi+itMi^^ where S,± are points 
just to the left and to the right of the arc respectively), then this average must vanish. The reason being, that this 
average represents the far-field felt by the charges on the arc. Looking under a magnifying glass at a segment of the 
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arc, one sees a long (from this perspective, infinite) chain of charges. These chains will fly off if an external (or far-) 
field is present. Namely, the average must vanish. These considerations allow finding /i(^). In fact Gaudin in a paper 
in Frenchi^ (later reviewed and expanded in English in [l3|) had shown that it is given by 

where 



R2K{i) = 



\ 



n(C-M,)^ + A,2 (8) 



Indeed, /i(^) defined by ([7]) has a jump continuity on the real axis of the given value p(e), vanishes on the endpoints 
of the arc has a non-trivial jump discontinuity on K arcs, and its average value across an arc is (since it simply 
changes sign across the arc). 

Eq. ([7]) is an expression for /i(^) given a knowledge of p{e) and of the endpoints of the arcs /Xj, , j = 1, . . . , if. /9(e) 
is arbitrary and may be tuned by blocking levels and placing rapidities between adjacent unblocked levels. The arc 
endpoints must be determined self-consistently, however. These self-consistency conditions can be derived by noting 
that h{£) as defined by (jS]) must have the following asymptotic behavior as ^ — oo: /i(^) — 5- ^ + O ^0 . Expanding in 

large ^, Eq. Q shows that the expected asymptotic behavior oih{^) is only satisfied if the following K —1 conditions 
hold: 

p{e)e'- 1 

^ -de = -<5/,if-i, l<K-l. (9) 



R2K{e) 9 

These are not enough to determine 2K free parameters which determine the position of the endpoints. Extra conditions 
may be found if one knows the number of rapidities on each one of the arcs. In the case of one arc, the number of 
rapidities on the arc is known if one knows the total number of electrons. Indeed, the total number of particles is 
2P -I- M, where P is the number of rapidities and M are the number of singly occupying electrons. The number of 
rapidities on the real axis is known since we know p{e) and so the number of rapidities on the arc is also known. If 
there is more than one arc, however, the total number of particles is not enough to fully determine the endpoints, and 
for the same number of particles, the same given p{£) and the same number of arcs, one may find different solutions 
depending on how many rapidities occupy each arc. The solutions differ by the location of the endpoints of the arcs. 

Since the total number of particles is a good quantum number, which also factors into finding the endpoints of the 
arc, it is useful to have an expression for this quantity. In fact we shall want to compute Jz = S^^—^^^^—^. This quantity 

is directly related to the total number of particles, Jz — S ( ^'2 ^ ) ) ' '^^ ^^'^^ features in the asymptotics of 

h{^) as ^ ^> 00. Indeed, as ^ — )• 00, h{^) i + Expanding Eq. ([7]) for large ^, one obtains: 



where [/('^)]^ denotes the positive Laurent series of /(^) when expanded around infinity. 

The total energy may also be found. Instead of the total energy, E, we compute a slightly different but directly 
related quantity, f , defined as follows: 



(11) 



then h[0 ^ i + ^ + f 



The eigen-states described by Eq. ^ generalize the eigen-states found by BCS. States directly corresponding to 
the BCS eigen-states may be recovered, but in addition to the states, one finds states which have no BCS counterpart. 
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To obtain the BCS eigen-states one must assume K = I. The expression for the total number of particles, the energy 
and the constraint, Eqs. pU)) . and respectively specialize to: 



These expressions are the familiar BCS expressions, if p{e) is identified as ^l^i-L^ where n(e) is the occupation 
number of excitations at energy e. Indeed, blocking a level or inserting an E on the real axis change the energy of 
the system, and thus may be viewed as excitations. A is the size of the spectral gap and fi is the chemical potential 
of the condensate, to be determined self consistently, given the total number of particles in the system. It is not 
clear however that the Richardson wave function, Eq. becomes the BCS state (or the number projected version 
thereof) in the thermodynamic limit. In particular, the factorized form of the BCS state does not appear naturally in 
Eq. Q. We shall see, however, that when if = 1 all correlation functions factorize, thus confirming BCS results. We 
shall also derive the results for K = 2, where, we show, the factorizability of correlation functions no longer holds. 



IV. SPHERICAL AND ELLIPTICAL CASES 

In the above, we have shown how to obtain the continuum solution of the Richardson equations and how to relate 
the solution to the good quantum numbers Jz and £. This was based on the electrostatic solution given by Gaudiiti^. 
We are interested in finding expectation values. The simplest expectation value to find is that of (S^). Where 

^f = 4+Q,++4_c,,_-i (14) 

By Hellman-Feynman, (S'f) = j-§f-, where £ is given by (|12|) . More complicated expectation values are not simply 
computable by invoking Hellman-Feynman, still objects similar to ^ will appear in such computations. More 

precisely, we have ^ = S ^—5 + J^i, ^f^) ; ^^'^ the object that will repeatedly appear in subsequent calculations will 
be We turn, then, to the computation of this object. 

It is easy to obtain information on ^f^, by considering the potential function (/'(^), given by: 

m = E - ^^■) + ^ E - K) (15) 

j V 

Then, a coarse-grained times the density of £^'s is given by the jump discontinuity of at ^ — E^, divided 

by 27ri(5. We compute then the potential in two cases: the case of one arc, and the case of two arcs. 

Note that the solution ([7]) derives its algebraic properties from the function R2k{£,) defined in ([5]). R2k{0 f^'Ct 
defines a spherical double-sheeted algebraic Riemann surface in the one arc case an an elliptical surface in the two 
arc case. We shall use basic notions of algebraic geometry, especially in the two-arc or elliptical case, where direct 
manipulation becomes cumbersome, and theorems of uniqueness of functions on Riemann surfaces become a much 
more powerful route to obtain the results. 

We are interested in finding the change in electrostatic potential S(j){^) when one changes the charge density on the 
real axis. It will be easier to compute first 6h{^)d£, and then integrate. We first compute Sh{^)d£^ when a unit charge 
is added at point e on the real axis, and then obtain a general result through a simple application of the superposition 
principle. Based on its definition, Eq. ([5]), the field h{£,)d£^ will have a pole at e, with residue 1. Similarly, there will 
be a pole at infinity. Indeed, § h{S^)d£^ — S{N — ■§), where the integration contour encircles infinity. In addition, the 
arc will move, namely A^ and pi will change. If we look at Eq. ([7]) we see that, no matter how Aj and pi change, 
/i(^) will retain the properties that it is an algebraic function defined on the Riemann surface of R2k{x) and that 
it changes sign going from the upper sheet to the lower sheet. This means that Sh{£^)d^ will also have a pole at the 
lower sheet at e and at 00 but with residues having reversed signs, respectively. We have (5/i(^)d^ as a meromorphic 
differential with given pole structure. Such a differential is unique on a spherical Riemann surface, and in fact is given 
by: 



Shim = , ' 1 + ^^^7^^ + ^^ rfC (16) 

V(C-^)2 + A2 I ^-e y 
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To obtain then Sh(^) for a generic disturbance of the charge density on the real axis, one must employ the superposition 
principle: 



Our goal is to obtain ^q^P ■ To obtain this we must assume that Scj^iS,) in (jl7l) corresponds to an infinitesimal change 
produced by moving £j. It is important to realize, though, that when one moves a single Si the rapidities present near 
it on the real axis will also move, thus contributing to 5p. The computation of 5p when one moves a single Si is then 
very difficult, since one has to know the exact configuration of rapidities near e^. This is not necessary, however, since 
we will be interested in coarse grained quantities. Indeed consider changing not a single Si but rather a number, A, 
of them centered around e. It is easy to see that to leading order in i5, as we shift this group of Ei's, the whole charge 
density, including the rapidities 'trapped' between the SiS will shift rigidly. This implies the following equation for 
the change of charge occupation 5p: 

5p{e') = 6p{e)6'{e' - e)5e x 5, (18) 
where ^ is the amount each one of the A e^'s around e was shifted. Plugging this into (IT71) we obtain: 

m) , + {e - p){^ - p) 

-^=^-^(^)i?.(e)i?.(C)(s-e)- (^9) 

which upon integration yields: 

We wish now to obtain a similar result when two arcs are present. This requires working with the elliptic Riemann 
surface described by Ri{C)- Being elliptic, this Riemann surface has the topology of a torus. We may use a rectangle 
with points on opposite sides identified (cyclic boundary conditions) as a model for the torus. In other words, we 
shall write the results first on the rectangle and then map them to the algebraic curve defined by i?4(0- Consider 
then a rectangle of sides oJi G and —iu)2 G K^. A function that maps this torus into the two sheeted algebraic 
curve defined by i?4(^), is given by the inverse Abel map, written in terms of Weierstrass's elliptic Zeta function: 

^(u) = C^{u-Uoo |wi,a;2) - C (w + Uoo |wi,i:j2) + c, (21) 

while the direct map is given by: 

This procedure is standard in the study of algebraic Riemann surface, being reviewed in any number of textbooks. 

Consider now Sh{£^)d^ for the two-arc case. Being completely general, the pole structure is the same as in the 
one-arc case, and just as before starting from ([7]) we may conclude that Sh{(^)d^ is an elliptic differential that changes 
sign as one changes sheets. We need one more ingredient to find Sh{^)d^, as these conditions are not sufficient to 
determine a differential on an elliptic Riemann surface, due to the existence of the holomorphic differential, du{^). 
An extra condition on is obtained by considering that the change performed must leave the number of pairs on each 
arc constant. Now, the number of pairs on an arc is proportional to f Sh(^)d^, where the integral is taken around 
the arc, which serves as the required additional condition. The image of the arc on the rectangle is a line extending 
from to UJl. So we must demand that the integral of dh{^)d^ be zero taken on a line from to ioi. The final result 
is then: 

6h{0d^ = (23) 
duiO I { C HO - u{e)) - C HO - Woo) - C HO + <0) + C HO + ^oo) + iMflz^^^KM^ Spis)de 

Indeed, it is easy to verify that the expression inside the brackets in the integrand has poles at u{e) and Uqo, which are 
just the images of e and oo, respectively; The residues are correct; The integrand is invariant upon switching sheets 
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(this is effected by taking u(^) — s- — but the integral is multiphed by du, which does change signs when one 
switches sheets, so that the whole expression has the right behavior upon switching sheets; It is also easy to verify, 
that the integral of this expression taken from to uji gives zero. 

Plugging ([18)) into ([T7)l amounts to taking a derivative with respect to u{e), which yields: 

Shi^d^ = Pie) [p(«(0 - u{e)) + p(«(0 + u{e)) + 4uj-\ (^^)] du{OMe) (24) 
It is a matter of performing an integral w.r.t u(^) in order to obtain 

= p{s) [C(u(0 - u{s)) + C(«(0 + uie)) + Aoj^MOC ( y)] Su{s) (25) 

Using the identity: 

[C(«(0 - u{e)) + C(«(0 + uie)) - ({uiO - u^) - CHO + ^oo)] du{e) = , (26) 

i?4(£)(C - e) 

The expression for may be written in a form which will prove much more convenient later: 

Sm^,,Pie)Rm M (27) 



fe ^4 (e) V C ~ £ 

where 



= m) ■ ^ ^ 

V. SLAVNOV'S FORMULA 

In order to compute the matrix elements we shall make use of Slavnov's formula^. Two states will have a non-zero 
overlap only if they have the same seniority, namely the same set of singly occupied levels ji with the spins pointing 
in the same direction. We thus suppress the notation of seniority and write simply \{Eij}^^^) for a Richardson state. 
Slavnov's formula^' as applied^ to the Richardson solution reads: 

{{'w^},.=i\Mt.=i) = ri 7 ^rf 1 rdetJ, (29) 

llbKai'^^b - Wa) [ia^biVb ~ Va) 

where {v^}^^i obey the Richardson equations, while {'Wu}!^=i do not necessarily satisfy the Richardson equations. 
The matrix J appearing in (|29p is given by: 

Va-Wb I {Va - ea){Wb - £q) ^ (Wq - Vc){Vb - Vc) ' 

When the set coincides with the set {w^}f^^i, Slavnov's formula gives the norm of the Richardson state. In 

this case the matrix J takes the form: 

A ) J2a )2 ~ 2 Y^c^a {v^-vj^ a = b 

(Va-Vb)'' ' 

This limit form, of J will appear frequently in the sequel. 

VI. COMPUTATION OF EXPECTATION VALUES: BASIC EXAMPLES 



The notations for the computation of a general correlation function become quite cumbersome. It is easier to 
start with two simple examples which demonstrate the principle of the computation before plunging into the general 
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scheme. This is undertaken in the next two subsections, respectively. 



A. Computation of {Sz{e)) 

Consider the computation of {Sz{s)). is given by ([TUl . What is meant by Sz{e) is a coarse grained version of 
this quantity, namely 

= ^ E (32) 

\s^-e\<AS 

To compute such an object we may first compute simply (Ni) and then perform a coarse graining and subtract a 
constant to obtain {Sz{s)). 

We shall want to represent (Ni) as an overlap between two Richardson states, in order to use Slavnov's formula to 
compute it. More explicitly, we are computing {{v^}^^i\Ni\{v,y}f^^-^^) , where Ni — . We can write: 



^.iKir^i) = E^-^iKWa)- (33) 

OL 

Considering that the operator Ni simply projects on the space of states that have i occupied by a Cooper pair and 
inspecting Eq. it is quite easy to understand how to derive (I33p . We shall not give a more explicit proof, but 
rather explain the different ingredients on the right hand side of. First note, that in order for the level i to be occupied 
with a Cooper pair, one of the the operators 6^ in Eq. (|4]) must have hit the level i. The sum over a in p3|) is a sum 
over all the possible such a's. The factor ^ is inherited directly from h^^. b\ in (|33p is responsible for filling up 
the level i. All the rest of the levels have a chance to be filled by except v — a. This explains the state \{vy}i,^a) 
appearing in ([33]). This heuristic explanation may be translated directly into a rigorous proof, an alternative is to use 
the language of the algebraic Bethe ansatz to obtain the result, as done in 0and[i3. 

We shall denote the state b\\{vjj}y^a) by U {ej). We have: 

m^ALl) = E U {£.}) (34) 

a 

The state |{t'i/}iy^a U {£»}) = b\\{v,j}i^^a) can be thought as a Richardson state, with a set of rapidities, w^, which do 
not satisfy Richardson's equations. This can be done because of the following relation: 

\{vu}y^a U {si}) = lim e\{vy}^^a U {si + e}), (35) 

where the state \{vu}v^a U {e^ + e}) is given by (j!]). 

Having written (iVj) as a sum over overlaps between Richardson states, we are ready to use Slavnov's formula to 
compute it. The results is: 

,^'-r^ v-^detA(?) 

(^^)=E^i^' (36) 

a 

where A is given in ([5T|) and ^(") is A with column a replaced by a column vector, y , this column vector being 
given by: 



More explicitly: 



K^'^ = , ^ ,2 - (37) 

[Vy ~ EiY 



^'^ \ A^^y otherwise 
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(?) 

By Cramer's rule the ratio of determinants can be computed as '^''I^'a = ([^"^^^'■'la)- ^'^ *° ^^^'^ 

invert A, we note that A has in fact a straightforward interpretation. Suppose that each one of the Wj^'s are subjected 
to an external field Sh,y. In order for them to remain in electrostatic equilibrium they must obey the equations : 

S V - ^ V ^— -- = SK. (39) 

The shift of the u,y's in linear order is a matrix multiplying the vector 5h with components 5h^. This matrix turns 
out to be tA^^. Indeed expanding l\'39\i one obtains: 



Ai^jSvj = Shi. (40) 



Suppose we shift Sj by an amount Ssj . This change can be represented as having v^, experiencing an external field of 
dh^ = f (^/_!^^)2 = ^SejV^^K Which implies: 



which gives: 



To compute , we resort to the results of section IIVI where we have computed the coarse grained change in the 
potential due to a shift of a group of e's on the real axis. The jump discontinuity at of ^^p- is equal to 2'KiS times 

an averaged times the density of w's. This is true if Vfj, is on the arc or on the real axis, except that there is 
an extra contribution from unblocked e. If we integrate over the jump discontinuity we obtain the sum on the right 
hand side of (|42|) with two caveats - (1) the result is not a derivative with respect to a single Si but a coarsed grained 
quantity and (2) the unblocked e add to the result. The two caveats amount to the fact that by integrating over the 
jump discontinuity of ^^^p- over the arcs and the real axis, we obtain Sz{e) of definition ([5^ . Since the integral over 
the jump discontinuity is simply given by a contour integral surrounding both the real axis and the arcs, we have: 



In the spherical (one arc) case, is given by (|20|) . thus 

Joo R2(e){^ - e) 2n R2[e) 

which is the known BCS result, derived in this context by Richardson^, based on more direct methods than the use of 
the Slavnov formula - methods which are nevertheless harder to generalize to more complicated expectation values. 
We extend the result by considering two arced configurations. In this case, ^^p- is given by ([27| . thus 

m'(^) is computed making use of (|2T|) . to yield: 

^'(0"' = pHO - ^*oo) - pHO + "oo). (46) 



The first integral in (pS)) is to be taken over a large circle encompassing the arcs and e while the second integral is 
taken over the circle's image under u(^). Performing the integration is a straightforward exercise in picking up the 
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respective poles, the final result being: 



^ = (.-M.)(.-M,) + ^-2p(2u.)-4.r-a^) 

i?4(£) 



B. Computation of (£*)S{£)) 



We now compute {{v^}f^^i\Sl,Si\{v^}f^^i), assuming both i and i* are unblocked levels. Consider then 5*1. 5^1 {wi/}^^^) • 
Sjt projects \{v^}f^^i) on the space in which Si- is empty and then fills it. The effect of Sj, can be simply achieved 
by adding Si* to the set of v^, as this causes the level i* to be filled while ensuring that no Vi, hits e^. . Namely, 

SJUMLi) = S^lMLi U {£..})• (48) 
Si projects U {si*}) on the space in which Si is filled and then empties it. In formulas: 

S^\MLl U {£..}) = bM\ML, U {£..}) = Y^^^lMLi \ M U {e^',e.}), (49) 



where in the last equality we have used the representation of Ni as an overlap, Eq. p4p . The following identity is 
easy to understand: 

hlMLi \ M u {£,p , £j) - (1 - n,)\MLi \ M u = (50) 

making use again in the last equality of the representation of Ni as an overlap. We obtain: 

(K}r.ii^.5i.iK}r=i>= (51) 



{Mu=l\M,,^l \ {Vo,,Vp} U , e^}) 



Making Use Slavnov's formula, we are now ready to write the expectation value {Sj^Si) as a determinant. The first 
term on the right hand side of (ISTt goes along the same lines as the computation of (Ni), so we shall not repeat it 
here. The second term on the right hand side of (|51|l has the added feature that it has two replacements Va £1' 
and vp — > £i. This leads to the following equation: 



{{vy}],=i\{vu}^^i\{va,vp}yj {£i.,ei\) ^ --^ detAVz* (52 

{Va-V0){£i- £i*) 



where 



A],:: = { (53) 



A^j_^i, otherwise 



Cramer's rule for A^i' i > reads 



detA - '^"^ I (A- v«)„ {A-^y'^^))^ ) ' ^^^^ 
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which according to pT|) reads 



detAV 



detA 



= det I f^' 



Combining (ISTI) . ([5^ and ([55j) we obtain: 



01 S^)^ 




(55) 



(56) 



We now need to take the continuum hmit of the expression. This is achieved by coarse graining the quantities Si and 



SJ,. Exphcitly, the coarse graining reads S^{e*) — J2\s- 
in (|56p has the following continuum limit: 



-s*\<AS 



Sj, and 5(e) = 2^ E 



^\ei-e*\<A5 



E Vg - i 



dVa 



de 



puis*) 



(r - e*) &(!>{£,*) 

{^* - e) de* 2ni5 ' 



Si. The first sum 



(57) 



where the integral encircles the arcs, but no part of the real axis (except the intersection of the arc with the real axis). 
Pu{£*) denotes the average occupation of unblocked levels at e* . This term appears because the unblocked i's (and 
only them) must be summed over in the coarse graining procedure. Indeed i and i* are assumed to be unblocked in 
([56| and if either one is blocked the correlation function is obviously zero. Another point to note is that since the 
integral is taken over F the contribution of ^p- is neglected for real Va- However this poses no difficulty, since only 
real near £j. are affected by a change of Si* , and their contribution to the sum is suppressed by a factor {v^ — Si*). 
Namely, this contribution does not survive in the continuum limit. 

Treating now the continuum limit of the second, double, sum in (|56p we obtain: 



(vg - ei*){vi3 - Ej.) / dvg dvp dv^ dvp 



y 



-2p,(e) 



dei> dsi 

^ e*) dcj^iC) dC 
{(* - e) de* 2niS 



dsi dsi 

iC^e*){C-s* 



(dm dm 



r 



i^-C*)ie-e*) V de de* 



dm)dm \ _d^dc_ f . 

de de* J 2TTiS2mS' ^ ' 



The double integral on the right hand side is the obvious continuum limit of the left hand side. The single integral 
takes into account the contribution of ■y's near e,;, which is neglected in the double integral, which is performed over 
r. This contribution is naturally proportional to p«(e), the average occupation of u's near e (again being related 
to the fact that the w's move rigidly with the e's). The contribution of w's near e^* is suppressed by the factor 
{va — ei*){vp — ei") and thus need not be taken. 



Since the average (coarse-grained) charge, p is defined as p 
limit in the following form: 



pu-2p„ 
2 



we obtain an expression in the continuum 



{S\e)S{e))^2p{e 



r (e 



e*)dm) dC 
- e) de* 2m5 



{i-e*){c-e*) ( dm) dm) _ dm)dm \ 

rJr (^-r)(e-e*) V de de* de de* J 2niS 27riS 

(59) 



We shall not perform the integrals explicitly, since the result is not very illuminating. We shall proceed instead to 
giving the general expression for the expectation value of any number of operators. 



VII. COMPUTATION OF EXPECTATION VALUES: GENERAL FORMULA 



We now compute a general expectation value featuring any fixed (not scaling with j) number of operators. The 
first thing to do is to write such an expectation value as an overlap of Richardson states. This is done either by 
invoking concepts related to the algebraic Bethe ansatz, or using repeatedly the tricks of subsections IVI Al and IVI 5] 
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The result is: 



(60) 



i-r 



k=l m.i<m2<---<mk i/i,...,i'„+„,+fc ni=l(''^i'! ^il) Yil^li^t^n+l ^im,) Yii^li^i^n+k + l ^jl) 

The factor (— )'° comes from a straightforward inclusion-exclusion principle, or alternatively from expanding the 
product rij'^ill ~ whose origin is the same as the origin of (1 — Ni) appearing in (|50)) . 



The overlaps appearing in (|60p can be easily computed using Slavnov's formula, with the result 



f f 1 1^2 

det Avfel '=2 



n,„<„(efe™ - ^fej nm<„(^^'>™ - ^^<^J detA 



/ 1/1 1/2 ... 1/s N 

in which '^i '=2 ... fea / ig defined as: 



7,(5 



with 



^vfci fe2 ::: ^1) ^ | 

" ^ otherwise 



A ^ is given an electrostatic interpretation just as above to yield: 

"^-i^(fe.) 



/ 1^1 1/2 ... t/s ' 

det^vfei fe2 ... k,. 



detA 



= det 



= det— ^ 



And the whole expectation value has the following continuum limit version: 



fc=o lePk 



n 



27r 



n 



n 



i < j 



where 

Pk = {/ C {1,2, . . . ,2n + m}|/ = {ji, j2, ■ ■ ■ ,jk,n + 1, n + 2, . . . , 2n + m} , where 1 < ji < j2 < 



(61) 



(62) 



(63) 



(64) 



(65) 



< jk < n} . 



(66) 



The difficult object to compute in (1551) is det^^^ Eq. (P7|) shows that, up to a multiphcation of rows an 

columns by constant factors, the matrix ^^^^ has a part which is a Cauchy matrix, the Cauchy matrix being given 
by: 

C^,, = (67) 

4.7 - Ski 
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Indeed 

det [C + G],^., (68) 

where 

Gij = 9{Cj). (69) 

We use this fact to write: 



det 



The determinant of the Cauchy matrix is known to be given by : 



jT pisk^Rii^ ^ ^ ) ^ ^ (s,j-yCr'Gi,] (70) 



det (C.,) = "<---^^^;,-<--^^^'-;"^'--^^f"^--^] (71) 



The inverse of the Cauchy matrix is also known, this is given by: 



Note however, that this is the inverse of C when it is understood that the indices run only over the set /, namely: 

V(z,j)G/', Y.CrlCij=6i,j. (73) 



lei 

Making use of an algebraic identity: 



one obtains that the matrix fC~^G] . . is diadic: 

If r and s are column vectors and is a diadic matrix formed out of them, F = rs*, then det(l + F) = 1+ r^s. This 
allows us to write: 

where the contour, C, on the right hand side; encircles all ^I'a. We let this contour be composed of two parts, a large 
circle traversed counterclockwise around infinity, and a contour, F, traversed clockwise around the arcs, ~ ioo ~ 
The integral over the large circle can be done immediately by expanding its radius to infinity. This integral can be 
easily seen to be equal to 1. We are thus left only with the integral over the contour F: 
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where the integral is taken counterclockwise around the arcs. Combining (|65|) . (iTOl) . (ffTj) and ([77)1 we obtain: 



n p(^^) 

2n-\-7ri 

n 

'i=2n+l 



n 2 



-1 + 7, 

2n 

n 

i— n+l 



Ri{ei){x - - Ei^n) ^niS 



Ri{Ci)ix - Si) 



i?4(£i)(a; - - Si) 27ri(5 



The result has a convenient almost-factorized form. Performing the integrals over ^'s directly we obtains our main 
result: 



(79) 



"2n+m 



n pi^'^) 

i=l 
2n 

n 



X Si 



ff(2;) 

2n-\-m 

n 

i=2n+l 



n 1 



(Ail + M2 - (.T + - Si) 



Ri{ei) 



Riiej) 

X ~ Si 



(/xi + /Z2 - (a: + £i))(a; - £») 
i?4(£») 



da; 
2^' 



where g{x) is given by (pS)) and -R4(^) is given by ([S]). The result appears formidable, but consists only of a single 
contour integral over known functions. The computation of this integral involves finding residues of the integrand, 
which is a mechanical task, easily performed by mathematical software, such that more explicit expressions can be 
derived for given n and m. 



In case one arc vanishes, e.g., A2 — >■ 0, the function g[^) can be shown to take the limit g{£,) 



In this case 



the integral over x in (j78p can be taken by shrinking the contour of integration to a point, ^2- This amounts to a 
substitution x ^2, and gives the BCS result: 



(^(£fcJ5(£fcJ . . . SiekJS\ek„^,) . . . 5^(£fc,J5,(£fc,„+ J . . . ,%{ek,„^J) 



(80) 



-^^^2(£^)' 

2—1 ^ ^ 



i—2n 



R2{ei) 



Note however that this way to obtain the BCS result is not general. The point /12 on the real axis is special, and 
has the property that the far-field (/i(/i2 + ^0+) + /i(/i2 — jO+) = 0) vanishes at this point. Not all solutions with one 
arc obey this constraint. The general way to obtain the BCS result is to rather take expression for ^^^^ as a 

starting point. This amounts to taking G = in It is easy then to proceed since the integral over g{x) does not 
show up and in fact the one-arc version of (|78p takes the simplified form: 



(^(£fcJ5(£fcJ . . . 5(£fcJ5t(efc„^J . . . 5t(£fc,J5,(£fc,„+J . . . 5,(£fc,„+„)) = 

R2(^i) 



n pi^^mh-i 

i=l i=l ^ 



r ^2(£i)(?i - £i) 27ri(5 



n 

i— n+l 



^2(6) 



(81) 



r -R2(£i)(^i - Si-n) 27ri5 



\ 2n+m , J. 



i?2(C. 



R2{e^m 



,) 2ni5) 



The integrals can be explicitly taken to give ([50]) . 



VIII. CONCLUSION 



We have shown how to compute correlation functions in the thermodynamic limit of the Richardson model. We 
gave explicit results for the case of one arc and two arcs. The one arc results converge with the BCS result, as 
expected. The correlation functions factorize into independent factors corresponding to each one of the operators in 
the correlation function. In the two arc case the factorization property disappears. Instead, the result, Eq. (|79l) . 
is given as a contour integral over an auxiliary variable x, which has a factorized form, where again each factor 
corresponds to an operator in the correlation function. 

It is interesting to see whether our results may also be obtained from a semi-classical approach, following the works 
of and[l6|. In this approach the semiclassical analoguesl^, 5+, S~ , Sz, of the operators , S, Sz, respectively, are 
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considered. These are shown to obey a classical integrable nonlinear equation. Being integrable, the equation may be 
solvedi^"— . It may then be possible to compute correlation functions in a semiclassical approach. In the case where the 
order parameter A{t) is time independent, this approach converges with the BCS approach, producing correct results. 
In case the order parameter A(i) is time dependent, it may be more delicate to justify the semiclassical approach for 
the computation of expectation values. The question of validity notwithstanding, our final result is suggestive of such 
an approach. Indeed, it may be that by a change of variables the integral over x turns into an integral over time, the 
periodicity of the integration contour over x related to the periodicity of a semiclassical solution, in which case our 
result may turn simply into a time average of the product of the respective semiclassical spin components, 5"*", S~ , 

A more challenging task, one we intend to pursue in future studies, is the calculation of matrix elements. For 
example {v\c^j ^Cj^a\w) or {i^|c^ _|_cj between two different eigen-states, \v), and \w). These matrix elements are 

important in predicting the dynamics of Richardson's state and in revealing its quantum coherence properties in 
different physical situations. Indeed {v\c^j ^Cj^a\w) is related to the transition rate between state \v) and \w) under 
a perturbation a^j,iy Such objects appear in the computation of the Fermi golden rule rate due to, e.g., phonon 
scattering. Note that The object {v\c^j „Cj,cy\w) ^ doesn't have a natural semiclassical counterpart, as the operators 

CT' Cj,<j, are not simply related to 5*+, , Sz- The other matrix element mentioned, namely {v\Cj _^clj appears 
naturally when one attempts to compute the tunneling of pairs into a superconductor in state \v). Such a computation 
appears in treating the Josephson effect or Andreev reflection. A Josephson effect setup is the obvious choice to 
measure the time dependent order parameter A(t) in an experiment. 
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with A. Nahum, B. Spivak, O. Agam and D. Orgad. 



J. Bardeen, L. N. Cooper, and J. R. Schrieffer, Physical Review 10 8, 1175 (1957) 
J. Bardeen, L. N. Cooper, and J. R. Schrieffer, Physical Review 106, 162 (1957)' 
R. W. Richardson, Physics Letters 3, 277 (1963), ~" ~ ~ 



E. Bettelheim, Europhysics Letters 90, 67002 (2010), arXiv: 1002.2289 [cond-mat.supr-conf, 
Y. M. Gal'Perin, V. I. Kozub, and B. Z. Spivak, Journal of Low Tem perature Physics 50, 185 (1983)[ 
N. A. Slavnov, Theoretical and Mathematical Physics 79, 502 (1989), 



H.-Q. Zhou, J. Links, R. H. McKenzie, and M. D. Gould, ,Phys. Rev. B 65, 060502 (2002)] [arXFy:cond-mat/0106390 
R. Richardson, Nucl. Phys. A 52, 221 (1964) 

' R. W. Richardson, Journal of Mathematical Physics 18, 1802 (1977)', 

' A. Faribault, P. Calabrese, and J. S. Caux, Phys. Rev. B 77, 064503 (2008)1 | arXiv:0710.4865 l 
L. Amico and A. Osterloh, Physical Review Letters 88, 127003 (2002), arXiv:cond-mat/0105141 
F. Braun and J. von Delft, Physical Review Letter s 81, 4712 (1998) , arXiv:cond-mat/9810146. 
' M. Gaudin, in Travaux de Michel Gaudin, Modeles Exactament Resolus (Les Editions de Physique, France, 1995) p. 247. 

J. M. Roman, G. Sierra, and J. Dukelsky, |Nuclear Ph ysics B 634, 48 3 (2002), arXiv:cond-mat/0202070', 
' E. A. Yuzbashyan, B. L. ^.Itshuler, V. B. Kuznetsov, and V. Z. Enolskii, 

[Journal of Physics A Mathematical General 38, 7831 (2005) , arXiv:cond-mat/0407501 



E. A. Yuzbashyan, B. L. Altshuler, V. B. Kuznetsov, and V. Z. Enolskii, Phys. Rev. B 72, 220503 (2005) 
|arXivxond-mat/ 0505493 , " 

R. A. Barankov, L. S. Levitov, and B. Z. Spivak, IPhysical Review Letters 93, 160401 (2004)[ [arXiv:cond-mat/0312053] 



